arXiv: 1505.07197v2 [math.AP] 26 Nov 2015 


Multi-layer asymptotic solution for wetting fronts in porous media 
with exponential moisture diffusivity 


By Christopher J. Budd and John M. Stockie 


We study the asymptotic behaviour of sharp front solutions arising from the nonlinear diffusion 
equation 6t = {D{9)6x)x, where the diffusivity is an exponential function D{6) = Doexp{f39). This 
problem arises for example in the study of unsaturated flow in porous media where 9 represents the 
liquid saturation. For physical parameters corresponding to actual porous media, the diffusivity at 
the residual saturation is D{B) = Do ^ 1 so that the diffusion problem is nearly degenerate. Such 
problems are characterised by wetting fronts that sharply delineate regions of saturated and unsat¬ 
urated flow, and that propagate with a well-defined speed. Using matched asymptotic expansions 
in the limit of large /3, we derive an analytical description of the solution that is uniformly valid 
throughout the wetting front. This is in contrast with most other related analyses that instead 
truncate the solution at some specific wetting front location, which is then calculated as part of the 
solution, and beyond that location the solution is undefined. Our asymptotic analysis demonstrates 
that the solution has a four-layer structure, and by matching through the adjacent layers we obtain 
an estimate of the wetting front location in terms of the material parameters describing the porous 
medium. Using numerical simulations of the original nonlinear diffusion equation, we demonstrate 
that the first few terms in our series solution provide approximations of physical quantities such 
as wetting front location and speed of propagation that are more accurate (over a wide range of 
admissible j3 values) than other asymptotic approximations reported in the literature. 


1. Introduction 


Problems with exponential diffusivity arise commonly in the study of water transport in variably- 
saturated porous media such as soil, rock or building materials [T]. An example of such, expressed 
in dimensionless form, is the one-dimensional nonlinear diffusion problem 



9{0,t) = 9i and 9{L,t) = 9o for t ^ 0, 
9{x, 0) = 9o for 0 < X < L, 


(la) 

(lb) 

(l c) 


where the solution 9{x,t) is called saturation and represents the volume fraction of pore space 
occupied by liquid. We are concerned here with the case when the dimensionless variable L S> 1 
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and the diffusivity is an exponential function of the solution having the form 

D{9) = Doe^\ (Id) 

where Do and /3 are constants satisfying 0 < Do <C 1 and /3 S> 1. This is a reasonable model 
for horizontal infiltration problems such as that pictured in Figured^, where liquid is taken up by 
capillary action within a dry porous medium and gravity can be neglected. Many experimental 
studies of porous media have been performed in which an exponential diffusion ansatz provides a 
good fit with measured data, mostly in the context of water transport in soil and rock [DIllllIl], 
but also for other porous materials such as wood [5], brick or concrete [?]■ Exponential diffusion 
has also been identified in other transport phenomena as diverse as heat conduction [8], optical 
lithography [9], solvent transport in polymers |10] . and diffusion of impurities in oxides m- The 
significance of this type of diffusion coefficient is that if 6i > 9o and if /3 is even moderately large 
then D{9) varies rapidly with 9, thereby causing solutions of (fTa)l to develop a steep interface, 
called a wetting front in the context of porous media flow, that is associated with localised high 
curvature and an exponential change in the solution gradient (see Figure [TJ)). The aim of this paper 
is to perform an asymptotic analysis of this phenomenon and in particular to study self-similar 
solutions of (fTa|l using a multi-layer asymptotic expansion, supported by numerical calculations. 
The asymptotic theory is especially subtle owing to the exponential change in the solution gradient 
and yields sharp estimates that agree well with numerical simulations. 



Figure 1. (a) A cylindrical porous medium with the left end immersed in a water reservoir, 
depicting the progress of a wetting front to the right owing to capillary action, (b) Plot of rescaled 
saturation 9 = {9 — 9o)/{9i — 9o) along the length of the cylinder from x = 0 to L. 


1.1. Overview of Previous Work 

The nonlinear diffusion equation (Hall has been studied in considerable detail for diffusivities D[9) 
having a variety of functional forms. Particular emphasis has been placed on the special case of 
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a power law, D{0) = DqO^ with m > 0, for which (|lap is called the porous medium equation 
or PME; an extensive literature exists for the PME that is thoroughly covered in the review by 
Vazquez m- In a typical wetting scenario the PME is supplemented by boundary conditions 
9{0,t) = 1 and 9{oo,t) = 0, in which case the solution for a power-law diffusivity is well known to 
have compact support and to consist of a front propagating to the right with speed proportional to 
Ahead of the front, the solution is identically zero and when m ^ 1 there is a discontinuity 
in the first derivative 9x at the point where the front meets the x-axis; otherwise, when 0 < m < 1, 
the wetting front meets the leading edge solution smoothly. The trailing edge of the front, on 
the other hand, always experiences a smooth transition similar to that shown in the saturation 
profile in Eigure [TJj. Analytical results have been derived for other types of diffusion coefficient, 
for example by Shampine m who provides an existence-uniqueness proof for a general class of 
diffusivity functions with D{9) required to be a continuously differentiable function on 0 ^ ^ 1, 

satisfying D{0) = 0 and D(9) > 0 when 0 > 0. This and other analytical studies are characterised 
by the fact that they pertain to degenerate diffusion for which the diffusivity vanishes identically 
at zero saturation. 

In contrast with this previous work, the diffusion coefficient we consider in this paper is an 
exponential function of saturation that, although small, is still bounded away from zero; hence the 
solution remains everywhere classical and all disturbances propagate with infinite speed, analogous 
to solutions of the “usual” heat equation. When the PME wetting scenario described above is 
repeated for an exponential diffusivity, the solution remains non-zero for all times t > 0, even when 
the initial conditions have bounded support. Nonetheless, the problem can be nearly degenerate in 
the sense that D{0) 1 whereas D{9) = 0(1) for 9 away from zero. Eor example, typical parameter 

values for a porous soil have been estimated as Do ~ 2x 10“® mf js and /3 ~ 20 [4]. Eor such diffusion 
coefficients, the solution inherits features that are qualitatively similar to the degenerate problem, 
most notably a steep wetting front that propagates with finite speed and that is associated with 
localised high curvature (see Eigure [21(b)). We define the wetting front location x*(f) to be the point 
where the curvature of 9(x,t) is greatest, and for x > x* we have that 9^ <C 1. One feature that 
distinguishes the exponential diffusion equation from the PME is the fact that 9 exhibits more rapid 
variation than for a power law, which in turn has a substantial impact on the shape of the solution 
close to the wetting front. A direct comparison is afforded by Leech et al.’s study of concrete [7j 
wherein they use experimental data to fit both power-law and exponential diffusivities, yielding 
respectively D(9) = Do9^ and D(9) = Dge^^. In Eigure [21 we plot the diffusion coefficients and 
corresponding saturation profiles for these two choices of D(9). While the qualitative features of 
both solutions are similar, there is a significant difference in both the steepness and location of the 
wetting front, as well as the sharpness of the corner (refer to the zoomed-in region in Eigure [2](b)). 
Another distinguishing feature is that the solution in the exponential case is characterized by a 
small nonzero saturation ahead of the wetting front due to the fact that D(0) 0. 

Eollowing up on these experimental studies, several authors have derived theoretical results for 
the exponential diffusion problem. Crank’s book m provides a comprehensive treatment of an¬ 
alytical solutions for nonlinear diffusion equations circa 1975 with many forms of the diffusivity 
function. In particular, Crank derives a similarity solution for the exponential diffusion case (fol¬ 
lowing the work of Cooper [8]) that reduces the problem to a second-order ordinary differential 
equation (ODE) which he then solves numerically - indeed, this is the same ODE that we will 
present later in Section [3l An alternate approach using functional iteration has been developed 
based on an integral formulation of the governing equations by Parslow et al. |15| . 

In contrast with these iterative or numerical solution methods, we are interested here in de¬ 
veloping an asymptotic series representation of the solution. One study of particular relevance is 
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(a) Rescaled diffusion coefRcients (b) Saturation profiles 




Figure 2. Comparison of solutions to the nonlinear diffusion problem for both power-law and 
exponential diffusion coefficients, D{9) = and D{9) = . (a) On the left is the rescaled 

diffusivity D{9)/Do. (b) On the right are the corresponding saturation profiles 9{x,t) computed 
numerically at some fixed time t > 0. The black arrow indicates the direction of the wetting front 
motion, from left to right. The green arrows within the inset plot denote the wetting front x* for 
each solution, located at the point of maximum curvature. 


due to Babu [16] who derived an asymptotic solution by making use of the simplifying assump¬ 
tion that both saturation and diffusivity drop to zero ahead of a certain wetting front location. 
Although Babu’s asymptotic estimates are reasonably accurate, we will demonstrate that his ap¬ 
proach of truncating the solution at the wetting front introduces significant errors that can be 
reduced by considering an alternate approach that incorporates the effect of the extremely small 
but still nonzero diffusivity values ahead of the front. Other alternate series expansions based on 
an integral form of the exponential diffusion equation were derived by Parlange and co-workers in 
imiiB]. Finally, we point out a connection to the work of Elliott et al. m who applied methods 
from singular perturbation theory to analyse the detailed structure of the transition region for the 
power-law diffusivity D{9) = 009"^ in the limiting case of m —>■ oo, which is known as the mesa 
problem. Although this problem is still degenerate, the diffusivity in the large-m limit experiences 
a rate of increase in 9 approaching that of an exponential function; therefore, the analysis for the 
mesa problem can be considered as a prototype for problems such as ©• Asymptotic solutions were 
developed for a semi-conductor dopant diffusion problem by King [20] for the two limits m —0 
(constant D) and m ^ oo (mesa), and King and Please [2T] applied singular perturbation theory 
to obtain a matched asymptotic solution for the case m = 1 (linear D)] however, the asymptotic 
structure of the solution for the exponential diffusion problem has not yet been explored in similar 
detail. 


1.2. Summary of Main Results 

We will demonstrate in this paper that a multi-layer asymptotic expansion is capable of yielding an 
accurate estimate of the wetting front location, provided that the exponentially varying solution is 
properly resolved close to the front. In particular, we will show that there is a self-similar solution 
Q{y) that can be written in terms of the dimensionless similarity variable 

X 

































Multi-layer asymptotic solution for wetting fronts in porous media with exponential moisture diffusivity 


5 


The corresponding wetting front location is a constant and has the asymptotic expansion 


The large parameter 7 
expression 


1 1 11 

7 27 ^ 127 ® 


+ 0 



—0'(O) is related to the physical constant (3 through the asymptotic 


(3 — (3 {9i — 9o) — ^ + O ^ 

where 0:3 ~ is obtained numerically. Our asymptotic solution is distinct from other approxima¬ 
tions derived in the literature in that it provides insight into the detailed structure of the wetting 
front, as well as yielding estimates of quantities such as y* that are accurate over a wide range of 
parameters. The estimate of wetting front location could be of particular interest to engineers since 
it expresses an easily measurable quantity in terms of physical parameters. 

The remainder of this paper is structured as follows. In Section [2] we motivate the problem 
under study by using the example of water transport in unsaturated porous media, although we 
stress that these results are equally applicable to other nonlinear diffusion problems having an 
exponential diffusivity. In Section [3] we introduce a similarity transformation that permits us to 
recast problem as an ODE initial value problem, for which 7 ^ 1 is related to the magnitude of the 
initial slope of the similarity solution. Section 0] derives our key results, consisting of a four-layer 
asymptotic expansion for the similarity solution, expressed as a series in 7 on each layer. Matching 
the asymptotic expressions then yields series approximations for the quantities of physical interest 
such as the saturation and curvature at the sharp corner in the wetting front, as well as the location 
x-f of the front itself. In Section [ 6 ] we perform a detailed comparison of our asymptotic solution 
to other approximations in the literature, as well as validating the results with careful numerical 
simulations of the original governing partial differential equation. 


2. Physical Background 


Water transport in a saturated porous medium is well-known to obey Darcy’s law |22] . U = 
which states simply that the liquid velocity U (in m/s) is proportional to the gradient of total 
hydraulic potential ^ (in m). The proportionality constant K is known as the hydraulic conduc¬ 
tivity and has units of m/s. However, many porous media flows are unsaturated, which means 
that the pore volume is only partially filled with liquid, and in this case it is necessary to intro¬ 
duce the liquid volume fraction or saturation, 9. The conductivity in a variably saturated porous 
medium is typically assumed to depend on the local saturation, leading to an extended Darcy’s law, 
U = —K{9)'S/^, that includes a saturation-dependence in the hydraulic conductivity. When this 
expression for velocity is substituted into the continuity equation 


we obtain an evolution equation for 9 


m 


+ V-U = 0, 


- = V-(X(0)V<h), 


( 2 ) 


which is known as the Richards equation. If both and K are assumed to be single-valued functions 
of saturation, then we can define D{9) := K{9) {d^/d9) after which Eq. (I2|) reduces to the familiar 
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nonlinear diffusion equation 

ao 

- = V • {D{e)V9), 

for which Eq. (|lal) is the ID version. The function D{9) is referred to as the moisture diffusivity 
and has units of m? js. As mentioned in the Introduction, the exponential form (jldjl for diffusivity 
is motivated by the study of certain soils and porous building materials, for which an exponential 
function provides a good ht to experimental data. 

We consider an idealised, cylindrical geometry depicted in Figure [1^ that is consistent with 
the samples typically used in experimental studies of moisture transport. The porous cylinder is 
initially dry and has one end placed in a water reservoir. We are interested in problems such as 
the horizontal infiltration scenario pictured, or where capillary forces dominate over gravity, so that 
water is absorbed into the porous medium by capillary action alone. Water progresses into the 
sample as a nearly planar wetting front, and so we can assume that the flow is uni-directional and 
governed by the one-dimensional diffusion equation (jlah . The length of the sample is denoted by 
L, which is presumed large in relation to the diameter so that L ^ 1. 

In order to treat a general class of wetting scenarios we take samples that are never completely 
dry, which corresponds to the usual situation wherein a porous medium undergoes re-wetting after an 
initial wetting/draining cycle. Consequently, there exists a non-zero minimum or residual saturation 
9 = 9o deriving from water that is trapped in isolated portions of the porous matrix and which 
cannot be displaced by capillary action. Furthermore, we introduce a maximum saturation 0* ^ 1 
that cannot be exceeded owing to micropores that are too small to allow water to penetrate, no 
matter how large the capillary force. Consequently, the saturation 9 satisfies Q<9o^9^9i^\, 
while the boundary and initial conditions are as specihed in (llbl) and (fTcl) . Water content is usually 
reported in the literature in terms of reduced saturation 

- 9-9o _ 9-9o 

9i -9o A0 ’ 

for which the boundary and initial conditions reduce to 

0(O,t) = l, 9{L,t) = 0 and 0(x,O)=O. (3) 

A picture of a typical wetting front is displayed in Figure [T )3 in terms of the rescaled saturation 
variable 9. 


3. Self-Similar Solution 

We are now interested in finding a self-similar solution of the nonlinear diffusion equation on the 
half-space x G [0, oo], when L —>■ oo, and which satisfies the Dirichlet boundary condition 0(0, t) = 9o 
constant. The problem ([T]) is invariant under the scaling transformation y = xt~^^‘^, which suggests 
seeking a solution of the form 9{x,t) = where the specific form of the similarity variable 


is chosen to eliminate certain constants from the ODE and boundary conditions. We then define 
the new dependent variable 

0 := or equivalently 0 := 


(5) 
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where 

P:=pAe = m-0o)- 

This form of solution is consistent with the second boundary condition m provided that 
L /^ 1, which motivates our taking L —)• oo in what follows. This simplifies the prob¬ 
lem further by eliminating the exponential diffusion coefficient (jld|) . leading to the following ODE 
boundary value problem for 0 (y): 

0©^' = —yQ' for 0 < y < oo, (6a) 

0(0) = 1, (6b) 

0(oo) = 0OO := e“^. (6c) 

For the porous media of interest here, /3 is significantly greater than one so that the parameter ©co 
satisfies 0 < ©oo <C 1 . 

Rather than solving Eqs. ([ 6 ]) directly it is convenient, for both the numerical and the asymptotic 
calculations, to reformulate this boundary value problem as an initial value problem, in which the 
boundary condition (| 6 c)) is replaced by a second initial condition of Neumann type given by 

0'(O) = /3(/>'(0) = - 7 , (IM) 

where —7 represents the initial slope which we will assume to be large. From this point on, we use 
dU) to refer to the initial value problem consisting of equations ([ 6 a|) . (I 6 bll and (| 6 cll) . The asymptotic 
analysis performed in this paper considers the limit of 7 —f 00 which we will see shortly is equivalent 
to taking ©oo —f 0. The quantity 7 is not known a priori in terms of the physical parameters, and 
so one of the primary results of this paper will be to first derive asymptotic expressions for ^ and 
y* in terms of 7 and to then invert them to give 7 in terms of the physical quantities. 

Before moving on, we address the suitability of applying alternate boundary conditions of Neu¬ 
mann (or flux) type. First of all, it follows from equations (I 6 ap and () 6 b(i . along with a later result 
indicating ©"( 00 ) ~ 0, that the Dirichlet problem above has an implicit no-flux (zero Neumann) 
condition at infinity. Hence, there is no need to consider this case separately. When the left 
boundary condition is replaced with a flux contition of the form 03 ,( 0 , t) = a, no similarity solution 
is permitted unless a = 0 , and in the zero-flux case there is only the trivial solution © = © 00 . 
Therefore, it is sufficient to focus on the Dirichlet problem in ([ 6 ]). 

3.1. Preliminary Numerical Simulations 

Before proceeding with the analysis, we present several plots that illustrate the asymptotic behaviour 
of the self-similar solution for large 7 , computed via numerical simulations of the initial value 
problem ([HI). We use a shooting algorithm wherein a value of j3 is chosen, and then 7 = —©'( 0 ) 
is updated iteratively until © is sufficiently close to the right hand boundary value ©00 = e~^. 
This approach is similar to that employed by others for the exponential diffusion problem [3 [23], 
and more details on our shooting algorithm are given in Section [6.11 For illustration purposes, we 
choose physical parameters corresponding to a typical soil water uptake experiment [3]. Taking 
limiting values of 60 = 0.04 and 0* = 0.43, and an unsaturated diffusivity j3 = 20.5, we obtain 
parameters A0 = 0.39, ^ = 8.0 and ©00 = exp(—8.0) = 3.4 x 10“'^, which clearly satisfies the 
requirement Qoo ^ 1- We remark that other experiments on water transport in a wide range of 
porous media (including soils, concrete, and other building materials) suggest that allowable values 
of /3 are restricted to a fairly narrow interval of 4 ^ /S ^ 9. The resulting numerical solution is 
depicted in the two rightmost plots in Figure [3Kc). The upper plot displays the computed similarity 
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solution 0 (y) while the lower plot shows the corresponding curves for reduced saturation 9{x,t) at 
ten equally-spaced time intervals, which are determined from 0 (y) by transforming back to physical 
variables Eq. &■ 

The computed solution for /3 = 8 exhibits a well-defined wetting front that manifests in the 6 plot 
as a narrow region with a steep slope located immediately behind a sharp corner. In the plot of 0 
on the other hand, the steep front has been eliminated by the exponential stretching transformation 
([5]), but the wetting front location is still identified with a sharp corner. A zoomed-in view of the 
wetting front is shown in the inset for /] = 8 and clearly indicates that the solution within the corner 
region transitions rapidly to a small value of saturation, but this transition remains smooth. 


(a) /3 = 2 (b) ^ = 4 (c) /3 = 8 

(7 = 1.166, 0OO = 0.135) (7 = 1.850, 0oo = 0.0183) (7 = 2.736, 0oo = 0.000335) 



y V y 





Figure 3. Saturation for values of /3 = /3A0 = 2 (left), 4 (middle) and 8 (right). The top row shows 
the similarity variable 0(y) obtained by solving Eqs. ([ 6 ]) numerically. The bottom row contains 
corresponding plots of the reduced saturation, 9{x,t) = 1 -|- (log0)/,5, with time t taken at ten 
equally-spaced points. The inset at the top right is a zoomed-in view of the sharp corner. 


Additional pairs of solution plots are given in Figures [3](a,left) and [3](b,middle) for values of 
/3 = 2 and 4 respectively. Clearly, taking j3 smaller (or equivalently 0oo larger) causes the wetting 
front to exhibit a more gradual slope and a milder transition through the corner region, and in 
the extreme case oi j3 = 2 there is hardly any evidence of a wetting front at all. However, as we 
have indicated above, most problems of physical interest correspond to values of ,5 ^ 4 and this 
observation has an important effect on the accuracy of our asymptotic solution derived in Section HI 
We next investigate in more detail the effect on the similarity solution 0(y) of changes in j3. In 
particular. Figure [3] indicates that as ^ increases the value of 7 likewise increases, which leads to a 
steeper initial slope and a corresponding shift of the wetting front location toward the origin along 
the y-axis. A number of additional simulations are performed for 7 lying in the interval [0.5, 5.0] and 
the corresponding values of (3 are plotted in Figure [D(a). By performing a least-squares polynomial 
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fit to the computed points, we obtain to a very good approximation the quadratic polynomial ht 

/3 = -log0oo ~ 

which when displayed in Figure 01^ a) alongside the computed data is nearly indistinguishable. This 
relationship between /3 and 7 will be verified later in Section 0] when it is derived as part of our 
asymptotic solution. 

(a) Initial slope, 7 = —0'(O) (b) Front location, ?/* 




Figure 4. Left: Values of /3 and 7 = —0'(O) obtained from numerical simulations of Eq. ([ 6 ]) with 
parameters chosen as in Figure [3j The quadratic fit /S = 7 ^ + 1/2 is shown as a solid line for 
comparison purposes. Right: The computed wetting front location (estimated using the point of 
maximum curvature) is depicted along with the two-term asymptotic approximation from Eq. (jl5p . 


It is evident from the plots in Eigure [3] that for j3 sufficiently large there exists a sharp corner 
in the similarity solution Q{y) that can be identified with the wetting front location in plots of 
saturation 8 . In contrast with the power-law diffusion problem, where the wetting front is identified 
with a discontinuity in the solution derivative, the exponential diffusion problem exhibits a smooth 
transition through the front and so there is no unique front position. We therefore choose to identify 
the wetting front location y* with the point of maximum curvature in 0 which satisfies Q"{y^,) = 0 . 
We provide a preview of our asymptotic results in Eigure 01(b), which depicts the computed front 
location y* as a function of (3, alongside our our two-term asymptotic approximation of y* derived 
later in Section [4.21 Clearly, the analytical results are quite accurate for 13 in the physical range. 


4. Multi-Layer Asymptotic Solution 
4-1- Overview 

In this section, we will derive the asymptotic form of the solution for large 7 . The preliminary 
numerical results already shown in Figure [3] suggest that the solution for large 13 (and hence large 
7 ) can be separated into four regions or layers: 

(a) An linear inner solution in which 0(y) is close to linear and has gradient close to I/ 7 . This 
region extends over a range of y values lying between 0 and slightly below I /7 where the 
solution is determined by the initial conditions at y = 0 . 
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(b) A nearly constant outer solution on the right of y*, where the solution is determined by the 
far field condition as y —)• oo. Here the solution approaches 0 oo for y sufficiently large (in fact, 
for y only a little greater than ?/*). 

(c) An intermediate-range solution containing the wetting front at y* and of width Ay = 0(1/7^). 
In this region we see the rapid transition through the corner of the wetting front and the 
solution has (locally) and exponentially (in 7) large curvature. In the intermediate range, 
we see a transition from the linear solution, which is of order 1/7^ when y = I/7, to one 
exponentially small in 7 for y > y^ = I/7 + 1/27^ + O (I/7®). It will prove convenient to 
divide this intermediate layer into a left-range and a right-range to capture this behaviour to 
high order. 

These regions are depicted schematically in Figure [ 5 ] in terms of the similarity variable 0 (y). We 
expect that our asymptotic approximation will be inaccurate for small values of jd when the inner 
solution is more curved behind the front, but will improve as j 3 increases and the inner solution 
becomes closer to linear. 

0 



Figure 5 . The similarity solution 0 (y) is separated into several regions: an inner region on the 
left consisting of a nearly linear solution n(r); an outer region on the right where the solution 
is nearly constant; and an intermediate region near the wetting front y = y* within which an 
intermediate-range solution of locally high curvature connects the inner and outer solutions. 

Before going into the details, we first summarise the main result for the wetting front location 

1 I 11 / I \ 

W)' 

By taking y* as the point of maximum curvature in saturation 0 (y), we may then approximate the 
values of both saturation and curvature at the wetting front by 

0(y*) pa 6000 and 0 "(y*) « 

T 


( 8 ) 
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We will show further that 


0 



1 b - log ( 7 ) 

272 ^ 74 



where b 


11 

12 



Finally, the physical parameter /3 can be expressed as 




1 


as 


log ©oo — 7 + —2~*~^ 

2 7 ^ 


T 


(9) 


( 10 ) 


where as is a constant that we can estimate numerically as as = 1/12. It follows immediately that 
for large ^ we have 


^ _ 1 ^- 1 /^ - (f - 5 ^) + o ■ 


( 11 ) 


4-2. Inner Solution 

On the interval 0 ^ y < y*, we know that the solution slope is initially —7 and so it is natural to 
introduce a new independent variable 


r = jy, 

and define the inner solution as u{r) = 0(y). Rescaling Eqs. ([ 6 j) yields 


„ r u 

UU =- TT-, 

T 

( 12 a) 

u( 0 ) = 1 , 

( 12 b) 

u'( 0 ) = - 1 . 

( 12 c) 


The 1 / 7 ^ factor on the right hand side of (I12ap suggests that if 7 is large then we should use an 
asymptotic expansion of the form 


/ N / N uiir) U2(r) 

u{r) = uo{r) H- ^ -\ - ^ + O 


T 


T 


/yi. 


(13a) 


After substituting this expression into the ODE and boundary conditions for u, we obtain the 
following sequence of initial value problems up to O ( 7 ”^): 


Uq Uq" = 0 , 

uq ui" = -ruo, 

// // / 
Uq U2 = —UiUi — rUi , 


uo(0)=l, uo'(0) = -l, 

ui( 0 ) = 0 , ni'( 0 ) = 0 , 

U 2 ( 0 ) = 0 , U 2 '( 0 )= 0 . 


These problems can be integrated successively to obtain 


uo = l-r, 

-ui = ^ - ^(1 - rf + (1 - r) log(l - r), 


17 


1 


^2 = -Al-r)- -(1 - r) + —(1 - r) + 2 - -r log(l - r), 


12 


12 


(13b) 

(13c) 

(13d) 


which after substitution into (jl3al) gives the inner solution to O (7 ^). The 0(1) and O (7 
solutions are depicted in Eigure [ 6 ] for values of /S = 2, 4 and 8 . Plots of the numerical solution of 
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Eqs. ([ 6 l), computed using the shooting algorithm described in Section [ 6 ] (and labeled “Exact”) have 
also been included for comparison purposes. 

We now investigate more carefully the validity of these asymptotic expressions. It is immediately 
clear that for asymptotic regularity we should have 


/ N Ui(r) 

uo{r) > — 

7 

Taking the leading order term on each side of this equation yields 

( 14 ) 

Einally, we discuss how the first two terms in the asymptotic expansion (jl3p may be used to give 
a hrst estimate the wetting front location y*. To this end, we set uq + = 0 and neglect terms 

that are quadratic and logarithmic in (1 — r) to obtain 


r 



or 


y* 


1 1 

7 27 ^ ’ 


(15) 


which yields the first two terms in the front location ([7]). We also note that setting r = 1 yields the 
leading order estimates 


n(l) 


1 


and 11 ^( 1 ) = — 1 ) 


(16) 


which will be useful in later scaling arguments. One of the primary aims of the more careful 
asymptotic matching performed in Section [4.3.31 is to derive a correction to Eq. (|15l) that gives more 
refined estimates of u(l) and r* to verify the above rough calculation. 


(a) /3 = 2, 7 = 1.166 (b) /3 = 4, 7 = 1.850 (c) /3 = 8, 7 = 2.736 



Figure 6. The inner solution u{r) for values of /3 = 2, 4 and 8, showing the leading order ap¬ 
proximation uq, first order correction uq + ui/ 7 ^, and numerical solution of the ODE initial value 
problem. The point where the approximation uq touches the y-axis (A) corresponds to the leading 
order estimate of the wetting front location, y* ~ I/ 7 . 


4 . 3 . Intermediate Layer 

4 . 3 . 1 . Preliminary Estimate of Leading Order Asymptotics and Solution Scalings. As r ^ 1 we 
start to see the rapid transition between the linear solution and the exponentially small solution. 
We note that a standard application of the maximum principle in the intermediate region shows 
that u{r) > 0 and u' < 0 for all r, so that u —)• Uoo = ©cxd and u' —>■ 0 as r —>■ 00 . This transition 
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occurs over a narrow layer, which we will see is of width O (1/7^). Before getting into the detailed 
analysis of intermediate solution, we begin with some simple calculations that aim to establish the 
leading order asymptotic form of the solution in intermediate layer close to the wetting front at r* 
and also to estimate the width of this layer. These calculations will guide our more sophisticated 
calculations given later. Using the scalings in the inner region, we have 

TU^ 

uu” =- 5 -, u( 0 ) = 1 , n'(l“) = - 1 , u{oo) = Uoo = ©oo- ( 17 ) 

T 

If we now divide both sides of the ODE by ru, then integrate and apply the boundary conditions, 
we obtain 

-log(uoo) = 7 ^ / - -dr. ( 18 ) 

Jo ^ 

In the case when 7 is large, our previous estimates of the inner solution, combined with the numerical 
calculations, imply that: 


(a) u' PS —1 if 0 ^ r < r*; (b) u' « 0 if r > r*; (c) u" « 0 if r / r*; and (d) u''{r^) S> 1. 


It follows from these observations that 

fOO „,// 


7 


poo M prz M 

/ — dr « 7^ / — dr K, 

Jo r r 


rp r* 


which can be combined with (| 18 l) and the inner solution estimate r* = 1 + 0(1/7^) to obtain to 
leading order 


-logUoo ~ 7 ^- 


( 19 ) 


While this estimate for u^o is not precise, we can use this simple calculation as a guide in the 
subsequent development of the asymptotic theory. 

For our next estimate, we make the assumption that as u” is exponentially large only close to 
r*, the leading order behaviour of the integral in (| 18 p l can be approximated by freezing the value 
of r = r* in the denominator of this integral. By integrating and applying the previous result, we 
obtain the leading order equation 


7^u' 



This equation can be integrated exactly when r > 1 to obtain 


r* (r - 1) = 7^Uoo 




where 


li(z) 



dt 

logt 


( 20 ) 


( 21 ) 


is the logarithmic integral function. In order to determine the leading order behaviour in the 
intermediate layer r ~ 1 and to match with the inner and outer solutions, we consider the limit 
u{r)/uoo —f 1 . From (fT 6 ]l . we have that u(l) ^ 1/27^ and Uoo ~ exp(—7^), so that u{l)/uoo 1 
for large 7. To calculate the solution thus requires estimates of the logarithmic integral li(z) in the 
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two limits z —>■ 1 and z —>■ oo, which are well-known to be 

11 ( 2 ;) = Fe -I- log logz + O (log z) 

1>W- —+ —+ o(^ 


log 2 ; log^ 2 ; 


as 2 ; 


as 2 ; ^ 00 , 


where Pg = 0.5772156649... is the Euler-Mascheroni constant. Using the expansion 
z, we obtain 


( 22 ) 

(23) 

for large 


., 2 „. 1 ; /^“(l) \ _ 2 


7 liooli 


Ur 


= 7 


u{l) 


log(u(l)) - log(n 

G 


O 


log^(lioo) 


1 


+ 0 


(24) 


If we now let u/u^o be close to one, so that we are very close to the wetting front, and then 
combine the various expressions above, we have 


Because r* = 1 -|- O ( 1 / 7 ^), this is equivalent to 


-k log log(u/u 

00)^ ‘ 


1 


u = Uoo exp exp ( -7 ( r - 1 - ^ - O 1 1 - Pf 


1 


(25) 


(26) 


If (as consistent with our earlier estimates) we identify the term 1 -|- ^ + 0 in this expression 
with the wetting front at r*, we then have 


u = Uoo exp exp ( - ^ ^ - Pg ). (27) 

V TUoo 7 

The structure of the intermediate layer is now clear. Because u'^ « is very large indeed, there 
is a very rapid transition from the inner solution where tt' ~ — 1 to the solution u = u^o^u' = 0 as 
we pass through the wetting front, over a range Ar ~ 7 ^ 6 “'*' . 


4-3.2. Left Intermediate Solution. Guided by the above, we now refine the asymptotic calculation 
in the intermediate layer, in particular for the range r < 1 -|- 1 / 27 ^ -|- 0 ( 1 / 7 '^), to get a more precise 
estimate of the location of the wetting front and of the value of Uqo- This calculation is rather 
technical, but allows us to obtain higher order asymptotic estimates which can then be compared 
with the numerical calculations. Within this range, we first determine an appropriate rescaling 
of the spatial variable r and dependent variable u by studying the form of the inner solution as 
r —>■ 1“. In this limit, the dominant terms in the inner expansion ()13ll are 


u{r) = 1 — r -|—K 
T 



(28) 


Consequently, we make introduce a new independent variable s that is 0(1) and is given by 


r = 



(29) 
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We note that taking s large and negative (for example 0 ( 7 )) allows us to match to the inner 
solution. Furthermore, the results of the previous section show that there is very rapid convergence 
to the outer solution, with u Uoo, for s > 1/2 + 0 ( 1 / 7 ^). We are motivated by considerations of 
matching to define a new intermediate-range saturation variable v{s) via 


so that 


u{r) = 

T 


u'il ) = 00 ) and uil ) — 7 : imply r;( 0 ) 

27^ 2 


We then have to leading order 


vv = —- 


1 + 


T 


(30) 


As s increases the intermediate-range solution will match rapidly to the outer solution so that 


v{s) ^ Voo := j'^Uoo, v'{s)^0 if s > s* = 1/2-h 0 ( 1 / 7 ^). (31) 

Integrating Eq. (1501) over the interval [ 5 , 5 +) yields 

7-^(log(.(8)) - log(»»)) = I (1 / ■ 

We will assume that the the range of s is restricted to 

s/ 7 ^ < 1, 


so that the denominator in the integrand may be expanded as a geometric series and 
7 “^ log(i;(s)) - log(uoo)) = ^ v” (^1-^ + 0 ^ ds. 

Using integration by parts and applying the far field condition f 0 for s/", we find that 
7-2 (log(u) - log(uoo)) = - 1 ;' (^1 - + ^ (:^) > 

which can be rearranged to obtain 

v' = 7"^ (log(woo) - log(u)) + ^ + ^ (^1^^ • (^2) 

We now carefully consider the solution of this equation for s < 1/2 and match it to the inner 
solution. This calculation is rather technical and involved, however it furnishes very precise infor¬ 
mation about the value of 0oo and the wetting front location. To do this we develop an asymptotic 
expression of the form 

u(s) = uo(s)-b -b -h ... when s < ^-b 0(1/7^). (33) 

7 ^ 7 ^ 2 

A subtle feature of this expansion follows from the coexistence of the logarithmic and polynomial 
expressions in both the intermediate layer expansion and the inner solution. This means that in 
order to capture the solution behaviour it is necessary that fi(s) has terms of both 0(1} and 
O (logy). However, provided that 7 is large, this means that we still have uq ^ ^ 1 / 7 ^ 3> ^ 2 / 7 ^ so 
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that the formal nature of the asymptotic series is preserved. We assume further, motivated by the 
previous analysis, that 


„(0) = i + al^ + 4+Oa,, 

2 7^ 7^ \7^ ' 


log Uoo = - 7 ^ + c log 7 + d + O 


(34) 

(35) 


and we will determine explicit values for the constants a, b, c, d. We next substitute the expressions 
(f33|) ~ (f35]) into (f32]) and solve the equations arising at each order. 

Considering first the leading terms of O (1), we find that 


which has solution 


Vo' = - 1 , uo( 0 ) = 


^o(s) = 2 “ 


1 


(36) 


Next, take the terms of O ( 1 / 7 ^, log 7 / 7 ^) for which we first need to expand logu from Eq. (I32p as 
a series in 7 


logu = log ( { ^ - s ) + vi/j"^ + 


= log ( i - s) + O (^7 ) , 


where we have used the fact that uq > 0 when s < ll2. The O ( 1 / 7 ^, log 7 / 7 ^) terms in the 
intermediate layer expression then become 

vi = c log 7 + - log Q - > 

which can be integrated to obtain 

ui(s) =cslog 7 + ^d- ^ s + olog 7 + 6 + Q - log Q “ (3^) 

We now match the intermediate range expansion v = uo + ui/y^ + O (1/7*^) to the inner solution 
when s < 0 and \s\ = O ( 7 ) ^ 1. Making use of the expansion 

(i - s) log (5 - 7 (1 - *) ° (7 

we can write the intermediate-range solution with terms ordered by size (for this range of s) as 


1 


1 /I 


72 \2 


log 7 , 1 


= -9 log(-'S) + C S- 5 - + -7 7 + X S + 


1 


/y 2 ^2 ry 2 2 

Next, rewrite the inner solution (I13p using v = 7 ^u(r) and 


'-y2 ^2 

a ^7 + — + —-(1 + log2)-|-O (. (38) 


1 1 


r = 1 -|- 


,,2 ■ 
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Then, for \s\ = O ( 7 ), we obtain the inner expansion in terms of the intermediate range variable as 




7 w l + ^ = -- s-gS log(-s) + 2 s 


1 


T 


, log7 

7^ 


+ 


17 


1 1 


+ 77o - 


127^ Y 2 


log 7 

72 


+ oi- 


By comparing this expression with that for the intermediate range in (1381) . we find that all terms 
match (both constant and logarithmic) provided the constants satisfy the identities 

1 


a = — 1 , 


11 1 , 

b = -log 2 , 

12 2 ^ ’ 


c = 2 , 


d = -- 


Therefore, the intermediate-range solution for s < 1/2 is given asymptotically by 


, . ,1 \ , , log 7 11 1 

Hs) = ( - - sj + ( 2 s - 1 )^^ + + 


T 


1-7io6(1-»|+o 


T 


which we can write as 
vis) = 


1 11 
2 + 1 ^“" 


^ ^ -21og(7) + log(l/2-s) I 
72 


T 


Substituting the constants into 

Similarly, ([351) implies that 


, we also find that to leading order 
1 




or 




7" 


1 


log Uoo = -7 + 2 log 7 - - -h O ^ 


1 


7 


so that 


= ^2g-7=-l/2+0(l/7=) 


or 


0 _ _ „- 7 ^- 1 / 2 + 0 ( 1 / 7 ^) 

Woo - ^00 - ^ 


(39) 


(40) 


(41) 


(42) 


4 . 3 . 3 . A More Precise Estimate of Front Location and Solution Curvature. The expression (HU]) , 
in which v{s) would appear to vanish when s = 1/2 -|- 11 / 127 ^ -|- O ( 1 / 7 ^), is strongly suggestive of 
a wetting front location given in terms of the inner variable by 

= l + ^ + + • (43) 

To further support this estimate, we return to the expression m for the overall behaviour of the 
solution in the intermediate layer, and continue the calculation given earlier using the more rehned 
values for u(l) and u^o given by (I4TD and (142]) . Expanding the li( 2 ;) function to include the z/ log^ z 
term, we obtain after some manipulation 


2 ,• l'u{l) 

'y '^00 

■ u„ 


1 + J_ f 11 

2 ^ 7 = 1 12 


M-lo: 


log(l/272)+72 + l/2 
1 11 1 „ / 1 
^ 1274 




+ — T + O 
474 


/-yL 


r*(r - 1 ) = 'j'^Uoo li 



- I'^Uoo 



Thus, as 
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we have 


(r - 1 ) = 


272 + 1274 + 47^ + C* J ^ 


1 + 




-7 1*00 li 


Ur, 


11 


. (u{r) 


2 Y 127^ V 7 ^/ V Uoo 

We deduce that if r* is now as given in (143 p . then with this refined value 


r - r* = - 7 ^Uoo li 


u{r) 


Ur 


(44) 


Applying the previous reasoning, this leads to exactly the expression (1271) , with r* now given by the 
rehned approximation (|43p . 

An alternative dehnition of the wetting front location is that it is the point r** where the 
curvature u” takes its maximum value, so that = 0. We now show that this is equivalent to 

the value r* given above, and estimate the curvature at this point. Differentiating the underlying 
differential equation yields 

W" + uV'=-ru'Vy" - l*77^ 

so that imposing = 0 yields at r = r** 


ru 


u = — 


1 + ' 


Assuming that ^ 1 and u"{r^^) ^ 1, we have to leading order 

7 1 1 

u (r**) =- 5 -. 

T 

However, based on the intermediate range equation derived earlier, we also have 


7 ^u' = — log 


u 


Ur, 


and so it follows immediately that the maximum value of u" arises when 


log 


u 


= 1 or u(r**) = euoo- 


(45) 


Substituting these various approximations into the underlying ODE we find that to leading order 
the maximum value of u" at r = r** is given by 

e7=-i/2 

(46) 


// 


1 


max w 


y^euoo 


T 


As expected, this is very large indeed for even moderately large 7 . This incidentally leads to severe 
problems in any numerical scheme that attempts to compute for large 7 . 

The wetting front location may now be estimated by substituting (14511 into (1441) to give 

r** - r* = -y^Uoo li(e). 

Because li(e) is of order one, we have that r** and r* are identical to all polynomial orders. 
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To complete this calculation, we now match to the outer solution. We recall from (j27p that the 
outer form of the intermediate range solution is given by 


u = Uoo exp exp 


{r - r^) 

I'^Uoo 


(47) 


We observe from this expression that u is very close to Uoo if is only slightly larger than r*. In this 
range, the term exp (—(r — r*) — Tg) is very small, so that we may approximate (|47p by 


U = Uoo 


1 + exp 


{r - r^) 


(48) 


Within the outer region, we have that u u^o as r ^ oo. To match to the solution outer region, 
we take this to correspond to those r values for which u is close to and look for a solution that 
is a small perturbation from a constant, namely 


u{r) = Uoo{l + g{r)), 


(49) 


where both \g\, \g'\ 1 for y sufficiently large. After substituting this expression into (| 6 ap . we 

obtain 

—2 / Hoo \ // 

19= - 0+9)9 , 

r 

which can be approximated for small l^l by 

Uoo9'' = -l~‘^rg'. 


This equation can be integrated once to obtain 


9 = —A exp 


211007 ^ j ’ 


where A > 0 is a constant. Close to r* we then have 


9 


/ 


—A exp 


r*{r-rO \ 

Uool"^ ) 


(50) 


(51) 


Because r* = 1 to leading order, this expression for g' exactly matches the derivative of the expres¬ 
sion ([451) loiig a-s 


A = 



This completes the outer solution matching to leading order. 


(52) 


4-5. Asymptotic Expansion for j 


We have so far expressed all asymptotic expansions in terms of the large parameter 7 . However, 7 
is not actually known a priori for a physical problem and so instead it is preferable to write 7 in 
terms of the known parameter jl. Taking the result from (142 p . we have that 


^ = 7^ + 7 + 


«3 

7^ 


+ 0 


T 


(53) 
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where we have introduced the coefficient 03 that is yet to be determined. Neglecting the O (7 
terms yields a quadratic equation for 7 ^ whose solution can be expressed for large /3 as 


7 = - (Y + 3^) + o 


(54) 


The high fidelity numerical calculations reported by Amodio et al. [23] provide convincing numerical 
evidence that 


as = 


1 

12 ’ 


although we have no formal calculation as yet to confirm this result. 

4 . 6 . Summary of Composite Asymptotic Solution 

Here we present a concise summary of the leading order asymptotic approximations for the inner, 
intermediate and outer solutions derived in Eqs. (fTHD . (15^ . (HSIl . (pUI) . (ISHIl and expressing 
each in terms of the original similarity variables y and 0 (y) as follows: 


Inner: 

(0 ^ y < 7"^) 


Left intermediate: 
( 7 -^ <y <y^) 

Right intermediate: 

{y > y*) 

Outer: 

{y 00 ) 


1 


0(y) = l-^y + 


^ - ^(1 - lyf + (1 - 72/) log(l - 7y) 




1 


0 ( 2 /) = -{y*- y) 7 + log 7 + log(y* - y) 


7 


0 ( 2 /) ~ 0OO + 0OO exp -Te - 


+ 0 


/yO 


y-y* 

7000 


0 '( 2 /) 


1 e 




7 0c 


exp 


y^ -yl 

' 2000 


(55a) 

(55b) 

(55c) 

(55d) 


In these above expressions, 0oo = e ^ and 7 can be written in terms of (3 as 


(55e) 


This last equation may also be used to replace 7 in the expansion of the front location (1331) to yield 

9 .« + (y + ^) (55t) 

where as = (numerically). 


5. Other Asymptotic Approximations 

In this section, we present two alternate asymptotic solutions to the exponential diffusion problem 
that are derived using other methods. 
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Babu [E] studied a modified version of the exponential diffusion problem in which he assumed that 
beyond some distance y ^ y* the saturation is equal to the residual value. In other words, he solves 
a modified boundary value problem that replaces our condition at infinity (l6c]l with 


©(y*) = ©c 


and 


dQ 

dy 


0 as y —f y* 


30 


Here, the front location y* is an unknown constant that must be determined as part of the solution, 
which explains the need for the extra limiting condition in (j6cl' 1 that corresponds physically to 
requiring a zero water flux at y*. By expanding the solution as a series in y = (1 — ©oo)//S, Babu 
obtained the following approximation for the wetting front location 


y? = 


1 (n 13 

1 -j- —77 I - “h — 

3 ' \ 90 8 


rf' + 


(56a) 


where we use a superscript “B” to distinguish Babu’s solution. Neglect exponentially small terms 
in ©oo = this expression may be rewritten in terms of /S as 


y 


B 

* 



(56b) 


which may then be compared directly to our asymptotic approximation in (|55fh . Although the lead¬ 
ing order terms in O are identical, the difference in coefficients at the next order generates 

a significant discrepancy in the front locations. Babu also derived two further approximations that 
he refers to as second- and third-order expansions 


©^’2(y) 

©^>3(y) 




^ 2 V 3 4 1 5 

-4o" 


(56c) 

(56d) 


both of which are defined on the interval 0 ^ y ^ y^. When y > yf, the saturation in both cases 
satisfies ©^’^ = ©^’^ = ©c^ (corresponding to 0 = 0). Babu’s solution is simpler than the one 
derived in this paper in that it is a polynomial expansion in integer powers of y and requires no 
multi-layer matching. However, this simplicity comes at the expense of reduced accuracy as well 
as a lack of information about the detailed structure of the wetting front. Another disadvantage of 
Babu’s approach is that none of his approximations for © is continuous at y = y^. 


5.2. Parlange’s Solution 

Many approaches for solving the nonlinear diffusion equation analytically begin with the Bruce- 
Klute equation, D{0) = —% Jq y(«) da, which can be derived from ([^ by treating y as a function 
of saturation 9 [25] . Although this equation was originally used to help interpret experimental data, 
it has subsequently been exploited by many authors to derive analytical solutions of the nonlinear 
diffusion equation, most notably by Philip [261 ISZj. A related iterative solution approach was 
developed by Parlange [281 [I7[, who applied a number of simplifications to obtain an approximate 
formula for y in terms of integrals of the diffusivity. For the specific case of an exponential diffusivity, 
Parlange m derived an implicit asymptotic representation for the saturation that can be expressed 
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in our similarity variables as 

y = (1 - 0^ + 0^ log e^/p) . (57a) 

In a similar manner to Babu, Parlange also imposed the condition that 0 = 0oo = for all 
y ^ y^. Equating terms in Eq. (j57aj) gives the result 

y:=[j3-'/^+j3-^/^){l-2e-f). (57b) 

Intriguingly, this estimate involves an exponential term that would correspond to a “beyond all 
orders” contribution to the results in this paper. This estimate behaves asymptotically in the limit 
of large ^ as 

yf ~^-1/2 + ^-3/2^ (58) 

from which it is clear that the front location matches only at leading order with our estimate (I55f|) 
and that of Babu from (j56bp . 

Another related solution was derived from the Bruce-Klute equation by Parlange et ah |18j . 
who extended an earlier method of Heaslet and Alksne [29] for the power-law diffusivity to the 
exponential case. This solution also has representation for 0 that involves the exponential integral 
function, Ei(z). Because evaluating 0 requires inverting Ei, their solution is much more complicated 
and expensive to compute than ours and so we have not considered it here. 


6. Comparison of Asymptotic and Numerical Results 

We now compare the various asymptotic solutions presented in the preceding sections. We also 
validate the asymptotic results using numerical simulations of the initial value problem in Eqs. (| 6 j) 
based on an algorithm that is described next. 

6.1. Solution Algorithm 

The 0DE (I 6 ap for the similarity variable Q{y) is solved numerically over an interval y € [0, M] 
using the Matlab initial value solver odel5s. Because both the problem and our asymptotic 
results correspond to the semi-infinite interval [0, oo], we must choose M sufficiently large that any 
error arising from truncating the right-hand boundary does not pollute the solution in the interior; 
on the other hand, 0 (y) tends quite rapidly toward 0 oo as y increases beyond the front, and so 
M does not need not be taken much larger than y*. In practice, we have found that choosing 
M = (twice the leading order estimate of the wetting front location) provides a reasonable 

compromise between efficiency and accuracy. 

For an actual wetting scenario, we know the asymptotic saturation Q^o (from (5)^ but the value 
of 7 in the second initial condition is not known a priori. We therefore build the initial value solver 
into a shooting type algorithm, for which we guess the value 7 = —0'(O), integrate the saturation 
variable to 0(M), and then compare to the target value 0oo- The value of 7 is then modified using 
the bisection method and this integration procedure is iterated until the relative error in the right 
boundary condition satisfies |0(M) — 0oo|/0oo < TOL, where TOL is a given tolerance. The leading 
term 7 ~ from the asymptotic formula (I55el) provides a sufficiently accurate initial guess to 
begin the iteration, and choosing ODE tolerances of ABSTDL = RELTOL = and 0-tolerance of 

TOL = 10“® yields a computed solution that for all intents and purposes can be treated as an “exact 
solution” of the original problem. 
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A significant difficulty with this algorithmic approach is that for even moderate values of 7 the 
solution curvature near the wetting front is on the order of exp( 7 ^), which can be extremely large. 
This large curvature causes problems with the initial value solver we are using, and effectively limits 
the maximum value of }3 to roughly 20, which restricts 7 ^ 4.4. With this restricted range of 7 , we 
can still test the validity of the asymptotic formulae derived here but we are not able to see their 
true accuracy at larger 7 , as the asymptotic series we have derived only converge at a polynomial 
rate as 7 increases. Consequently, we will also report some results using a much more sophisticated 
numerical approach that is the subject of a related paper m and which employs a high-order 
Taylor expansion based boundary value solver coupled with mesh adaptivity. This method permits 
calculations up to 7 = 18, corresponding to P ^ 325 and 0oo ~ 1-18 x 

6.2. Saturation Profiles 

Plots of the multi-layer asymptotic solution determined by Eqs. (j55p are depicted in Figure [3 for 
values of /3 = 4, 8 and 16 alongside the corresponding plots of the computed solution using the 
shooting method described above (which can essentially be considered as an “exact solution”). The 
inner and intermediate-range solution are both displayed, and the loss of asymptotic validity of the 
inner solution due to the logarithmic term is evident as y ^ y*. We have included in all plots the 
corresponding asymptotic solutions of Babu |16] and Parlange m, and in each case a second plot of 
all curves in terms of the rescaled saturation 6 is given (lower plots), which accentuates differences in 
the wetting front location. Our asymptotic solution is clearly an improvement over Babu’s solution 
for all values of which we attribute in large part to the fact that Babu’s approximation truncates 
the saturation at some approximate wetting front location and ignores the details of the solution 
structure within and to the right of the front. 

As described above, the relatively slow convergence of our asymptotic solution, which is poly¬ 
nomial in 7 , and hence also in (3, means that our asymptotic approximation is not particularly 
accurate for smaller values of fi such as fi = 4. This is reflected in the mismatch between the 
left and right intermediate-range solutions at the wetting front location (represented by the blue 
triangular point). However, our asymptotic approximation improves significantly in accuracy as fi 
increases, and when we take fi = 16 the comparison between the computed and asymptotic results 
are good. The differences at even larger are difficult to visualize using saturation plots, and so 
we compare the solutions further in the next sections in terms of the estimates for wetting front 
location. 


6.3. Wetting Front Location, y* 

We next focus on calculations of the wetting front location y*, which we have defined in our deriva¬ 
tion to the be point where Q" is a maximum. A visual comparison is provided in Figure [ 8 ] in terms 
of the computed front location as a function of showing our asymptotic estimate from (j55f() . 
Babu’s estimate yf given in (j56d|) . and Parlange’s estimate y^ given in (f58]l . We also included 
numerical estimates of front location using the Matlab solver ode 15s, for which we were able to 
compute up to a maximum of /S ~ 20 before the ODE solver failed. The plot of error in front loca¬ 
tion (relative to the computed solution) clearly shows that although Parlange’s estimate is better 
than our asymptotic solution for some values of ^ ^ 10 , our result is consistently superior for larger 
/S; indeed, even our two-term asymptotic solution surpasses Parlange’s result when ^ 13. 

Although our simple numerical approach fails for > 20, a more sophisticated numerical al¬ 
gorithm has been implemented by Amodio et al. [Mj that yields accurate solutions for values of 
much larger. Their results for a much wider range of /3 are summarized in Table [H along with 
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(a) ^ = 4 
7 = 1.850 






Figure 7. Comparison of our asymptotic solution ( - ) to those of Babu ( - ) and Par- 

lange ( - ) for /S = 4, 8 and 16, displayed in terms of the similarity variable (left) and rescaled 

saturation (right). The corresponding approximations of the wetting front location are denoted by 
points lying on the y-axis for our asymptotics (A), Babu (O) and Parlange (□). 
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(a) Wetting front location (b) Error in y^, 




Figure 8 . Comparison of the various asymptotic estimates for wetting front location to the 
computed (“exact”) solution. 


the corresponding asymptotic results for yT^, y^ and y^. The superior accuracy of our asymptotic 
solution for large ^ when compared with Babu’s and Parlange’s approximations is evident upon 
comparing the computed value of y* to the three asymptotic approximations. 

Table 1 

Comparison of our asymptotic results to computations using the method of Amodio et al. [24j . 


7 


Computed results 
doo 

y* 

Asymptotic results 

asy B P 

y* y* y* 

2 

4.559435 

1.046797e-2 

0.571747 

0.54536 

0.51539 

0.57103 

6 

36.50238 

1.403505e-16 

0.16911 

0.1689165 

0.16759 

0.170050 

10 

100.5008 

2.254440e-44 

0.1005094 

0.1004949 

0.100205 

0.100743 

18 

324.5002 

1.178490e-141 

0.0556417 

0.0556410 

0.0055591 

0.055683 


6.4- Initial Slope, —7 

Finally, we investigate the accuracy of our asymptotic approximation for 7 in Eq. (I55el) . where —7 
is the initial slope in the similarity variable 0. No corresponding estimate is available from either 
Babu’s or Parlange’s solutions. For a range of (3, Figure [9] compares the value of 7 determined 
from shooting simulations with that obtained from the 1-, 2- and 3-term asymptotic estimates. The 
accuracy of our series approximation is evident for increasing values of (3, which is essential because 
7 is a required input for our numerical simulations and yet it is the parameter (3 (not 7 ) that is 
known a priori for any given physical wetting scenario. 

We finish by restating some of the numerical results obtained in [23] for values of 7 as high as 18 
(corresponding to /3 up to 324). Beyond this range, even this more sensitive algorithm was unable 
to converge owing to the size of 9” at the wetting front. The primary calculations of y* and 0oo 
over a range of values of 7 were then enhanced by using Richardson extrapolation. The primary 
contribution of [ 23 j was to provide clear numerical evidence to support the validity of the following 
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(a) Asymptotic approximations for 7 (b) Error in 7 approximations 




Figure 9. Comparison of the various asymptotic approximations of 7 based on retaining one, two 
or three terms from the series in Eq. (|55el) . The left hand plot compares the three asymptotic 
results to the value of 7 obtained using the shooting algorithm. The right hand plot depicts the 
corresponding absolute error in the asymptotic approximations, treating the shooting results as the 
“exact” solution. 


expressions 


and 


1 1 11 2.96 

7 2 ' y ^ 127® 7'^ ' 




log((9oo) ^ 7^ 


1 

+ 2 + 


1 0.089 

1272 74 ’ 


which extends the asymptotics derived in this paper by one additional term. These results are in 
full agreement with our asymptotics and also hint at a more refined asymptotic calculation. 


7. Conclusions 

In this paper we have performed a multi-layer asymptotic analysis of a nonlinear diffusion equation 
where the diffusion coefficient is an exponential function, D{6) = . We focus on an application 

to the study of wetting front formation in unsaturated porous media flow, although exponential 
diffusion also arises in a number of other problems in heat transport, optical lithography, polymer 
diffusion, etc. The original boundary value problem for liquid saturation is reformulated as an 
initial value problem which, although not strictly necessary (c.f. Babu’s solution m in Section [5T]1 . 
nonetheless permits an easy comparison to numerics using an initial value solver. Motivated by the 
structure of the wetting front for large values of /3, we used the method of matched asymptotics 
to derive a four-layer solution consisting of a different asymptotic series for each of four regions 
corresponding to the wetting, and the leading and trailing edges. Other previous approaches to 
deriving (approximate) analytical solutions have assumed that the reduced saturation is identically 
equal to zero at a finite wetting front location, which essentially ignores the structure of the sharp 
corner that appears at the front for large /3. In contrast, our asymptotic solution uncovers the 
detailed structure of the transition within the wetting front, where the solution has very high 
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curvature, and is considerably more accurate that other previous approaches for large /?. The 
asymptotic solution converges polynomially in /? and maintains a high degree of accuracy over a wide 
range of /3 values corresponding to physical porous media such as soils and rock, unlike many other 
approaches which are either more limited in their applicability or exhibit a logarithmic singularity 
near the wetting front as {3 becomes large. Our approach has the additional advantage that it 
expresses the solution in terms of explicit formulas instead of requiring numerical approximation 
of a differential equation (which is required in the approaches of Wagner [9] or Shampine [30]) or 
inverting the exponential integral function (as in Parlange et al.’s asymptotic solution m)- In 
particular, our estimate for y* may be of value to groundwater researchers and experimentalists 
since it provides an explicit formula for a measurable quantity (wetting front location) in terms of 
physical parameter values. 
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